Walk-on-Hemispheres first-passage algorithm

Due to the isomorphism between an electrostatic problem and the corresponding Brownian diffusion one, the induced charge density on a conducting surface by a charge is isomorphic to the first-passage probability of the diffusion initiated at the location of the charge. Based on the isomorphism, many diffusion algorithms such as “Walk-on-Spheres” (WOS), “Walk-on-Planes” and so on have been developed. Among them, for fast diffusion simulations WOS algorithm is generally applied with an \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document}ε-layer, which is used for diffusion convergence on the boundary but induces another error from the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document}ε-layer in addition to the intrinsic Monte Carlo error. However, for a finite flat boundary it is possible to terminate a diffusion process via “Walk-on-Hemispheres” (WOH) algorithm without the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document}ε-layer. In this paper, we implement and demonstrate this algorithm for the induced charge density distribution on parallel infinite planes when a unit charge is between the plates. In addition, we apply it to the mutual capacitance of two circular parallel plates. In both simulations, WOH algorithm shows much better performance than the previous WOS algorithm.

www.nature.com/scientificreports/ Here, is the Laplacian operator. Then the normal derivative of G( r 0 , r) on ∂� creates the harmonic measure 18 and any harmonic function u( r) in satisfies the boundary integral equation; Here, n is the normal vector inwards the domain. For the required Green's function G( r 0 , r) , the linear combination of electric potentials can be used. To invoke the axial symmetry, we put a negative unit charge and its image charges on z-axis as Fig. 1 so that the potential vanishes on the conducting surface. Setting the radius of hemisphere R and the height of the source charge from the origin d, we obtain the Green's function in the spherical coordinates (r, θ, φ) (here, θ from 0 to π is the angle from the positive z axis downward) given by Taking the partial derivative of G(r, θ , φ) with respect to r, we get at R Now, let γ be the ratio d/R and the boundary ∂� = ∂X ∪ ∂Y ( ∂X : hemisphere, ∂Y : disk) and we integrate over the hemispherical part. The cumulative induced charge density is obtained to be (1) �G( � r 0 , � r) = −δ( � r 0 − � r), when � r ∈ � G( � r 0 , � r) = 0, when � r ∈ ∂�.
(2) u( � r 0 ) = u(� r) ∂G(� r, � r 0 ) ∂n dS ∂� . (3) G(r, θ, φ) = 1 4π ∂G(r, θ, φ) ∂r Figure 1. Schematic diagram of WOH algorithm; the radius of the hemisphere R and a charge −q at z = d and its three image charges at z = −d , z = −R 2 /d , and z = R 2 /d to make the potential zero on the boundary. www.nature.com/scientificreports/ and the total induced charge on the hemisphere becomes By Eqs. (5) and (6), the conditional cumulative distribution with respect to the azimuthal angle θ of the first passage location on the hemisphere is given by Then the inverse transform of the cumulative distribution function is given by The formula (8) gives the exact sampling on the spherical part of the hemisphere. For sampling on the disk part, because of the complexity of the inverse transformation sampling, we use conformal map 15,16 to exchange the location of the disk and spherical part as shown in Figure 2.
Let the distance from the origin O to the charge q and q ′ be r and r ′ , respectively.
Then, the relation between r and r ′ should satisfy the Eq. (9). Thus, the position of the charge q ′ is specified as r ′ = 2R/(1 + γ ) . In addition, for convenience of calculation, a variable γ = 2r/R − 1 and γ ′ = R/r − 1 are introduced respectively to represent the distance ratio from the center of the disk and hemisphere.
In Figure 2, the potential of the transferred charge q ′ is defined by Eq. (10).
By introducing the azimuthal angle θ/2 of the inversion sphere by Eq. (11), the induced charge density on the disk surface σ (θ) can be written in terms of σ ′ (θ).
The relation of the cumulative charges of the disk and the spherical part is given by the following; www.nature.com/scientificreports/ For practical use, we first determine which surface the diffusion passes through between the spherical and the disk parts by Eq. (6), the total induced charge on the spherical part. If the spherical part is chosen, then we just sample the angle θ by Eq. (8). Otherwise, we compute γ ′ as shown in Fig. 2 and then apply Eq. (8) to sample the angle θ/2. To follow the actual likelihood in Eq. (12), we use the acceptance-rejection method [8] with acceptance probability sec(θ/2)/ √

2.
Induced charge distribution on parallel infinite planes. When a charge is located between two infinite parallel conductors, the analytic solution for the induced charge density is known as a series solution only [19][20][21] . For the sample of the corresponding diffusion in diffusion Monte Carlo simulations in this geometry, "Walk-on-Spheres" (WOS) algorithm and recently developed "Infinite Parallel Plates" (IPP) [9] algorithm are available. However, in WOS algorithm the diffusion sample on the parallel boundaries can be biased by the ε -layer if the layer is not thin enough to suppress the error from the layer, which is needed for convergence of the diffusion simulation 22 . In the other IPP algorithm case, we have to use a tabulation and compute the additional terms of series solution whenever the sampling is too close to the rejection criteria so that the algorithm is somewhat complicated 19 .
In this section, the sampling of a diffusion in the parallel planes boundary is performed by WOH algorithm. The initial position of the diffusion is located at the middle of two infinite planes as shown in Fig. 3. For the WOH diffusion, the radius is fixed to D and the direction of the disk boundary is toward the plane of the minimum distance. The first-passage distribution of the diffusion simulation is compared to the corresponding electrostatic analytic series solution given below 21 ; The result in Fig. 4 verifies that WOH algorithm provides the correct induced charge distribution. In order to obtain the above result, we performed 100 independent runs of 10 9 Monte Carlo (MC) steps, that is, via the total number of 10 11 simulated diffusion quasiparticles. The convergence of the Monte Carlo errors is given in Fig. 5. The linear regression has its slope of − 0.49739 with the correlation coefficient of − 0.99986 . It is noted that all the logarithms used in this paper are the decimal logarithm.
In WOS algorithm, for large enough MC steps the Monte Carlo error convergence does not exhibit the linearity due to the error from the ε-layer. If MC steps are large enough, the error from the layer becomes dominant 23,24 . To see the ε-layer error in Fig. 6, we perform the same simulation replacing WOH with the WOS algorithm with various ε-layers. Figure 6 clearly states that the error convergence is hindered by the ε-layer. For the case www.nature.com/scientificreports/ of ε = 10 −6 , it seems that the simulation result is not much affected by the layer and the MC intrinsic error is dominant. Although the ε-layer error can be reduced by making the layer smaller, it causes the logarithmic increase of the simulation time [23][24][25] . In addition, in Fig. 7 with the induced charge density on the parallel infinite planes, we investigate the runtimes of the two algorithms, WOH and WOS. The runtime of WOH algorithm is inserted as the blue dotted guideline. The runtimes of WOS algorithm are obtained with ε from 10 −6 to 10 −2 . The linear regression of the WOS algorithm runtimes has its slope around − 10.51 with the correlation coefficient of − 0.99996 . In the parallel infinite plane simulation, it is clearly verified that the WOH algorithm is much more efficient than the WOS algorithm with any practical choice of ε-layer.
Mutual capacitance of two parallel circular plates. In this section, we demonstrate WOH algorithm for the mutual capacitance of two parallel circular plates with various radii and compare the performances with the WOS algorithm.
The mutual capacitances C ij of multi-conductors are represented in the matrix form which shows the relations between the total charge Q i on the i-th conductor and the voltages V j on the j-th conductor among N conductors;  www.nature.com/scientificreports/ Here in this paper, as in Fig. 8 we have two conductors only and so N = 2 . We compute the total charge on the i-th conductor by integrating the surface charge density on it via the last passage algorithm 26,27 , which can compute the charge density at a specific point on a conductor. The surface charge at r 0 is given in terms of the last passage Green's function, g( r, r 0 ) , of radius a of the last-passage hemisphere and the probability, P(� r → ∞) , of going to infinity of the diffusion which is initiated at r 0 .
Here, the surface integration is over the last-passage hemisphere, ∂�. When the i-th conductor C i has unit voltage and the others are grounded, the surface charge density σ ii on C i and σ ij on C j are represented by;  www.nature.com/scientificreports/ where the probabilities, P(� r → C i ) and P(� r → C j ) , are of going to capacitor i and j of the diffusions which are initiated at r 0 respectively. For the two circular plate conductor, the elements of capacitance matrix C 11 and C 12 are given by integrating the corresponding charge density over the conducting surfaces [10].
Let the distance between the two capacitors D and their radius R as shown in Fig. 8. We compute the mutual capacitance for various geometries by changing the ratio R/D from 0.1 to 10 3 . The simulation results are shown in the Table 1 from 100 independent runs of 10 10 Monte Carlo steps. The analytic solutions in Table I are taken  from the references 27,28 .
The runtime comparison of WOH algorithm with WOS one for the mutual capacitances of the two parallel circular plates from 100 independent simulations of 10 7 MC steps is given in Table 2 with respect to the geometries given by the ratios of the radius to the separation R/D from 0.1 to 10 3 . With the WOS algorithm, we used the ε-layers, 10 −2 , 10 −3 , 10 −4 , 10 −5 and 10 −6 . The runtimes of WOS are not significantly changed with respect to the geometry. With WOH algorithm, the runtime increases as the R/D ratio decreases and is less than WOS algorithm except the case of R/D = 0.1 and ε = 10 −2 .
The runtime data are obtained from 100 independent runs of 10 1 0 Monte Carlo steps. All computations were performed on a MPI PC cluster (15 nodes, 160 cores with 2.40 GHz and 120 cores with 2.10 GHz) with scalable parallel random number generator (SPRNG) 29 . (17)

Discussion
Due to the isomorphism between a Brownian diffusion problem and the corresponding electrostatic one, the induced charge density on a conducting surface by a charge inside the boundary can be obtained by the firstpassage probability of the diffusion initiated from the charge location and vice versa. In this isomorphism, an absorbing boundary surface in the diffusion problem corresponds to the conducting surface in the electrostatic one, a diffusion-starting position to the charge location, first-passage distribution on the absorbing boundary to the induced charge distribution on the conducting surface respectively. Based on the isomorphism, fast diffusion Monte Carlo algorithms have been developed, such as "Walk-on-Spheres" (WOS) algorithm, "Walk-on-Planes" (WOP) algorithm, "Walk-on-Hemispheres" (WOH) algorithm, "Walk-on-Cubes" algorithm and so on 3,5-8,10,30 . Among them, the simplest WOS algorithm is generally used with an ε-layer due to the fact that WOS algorithm can be used in any geometrical boundaries. However, it should be noted that the ε-layer for the convergence of diffusion induces an error 23,24 . Without the ε-layer, for an infinite flat boundary it is possible to terminate the diffusion process via WOP algorithm 7,8 . For a finite flat boundary or in two-plate between boundaries, WOP algorithm can not be used and WOH one is very plausible. In the previous works 4,12,13 , WOH algorithm was implemented via Von Neumann's accecptance-rejection method 14 . In this paper, via a conformal map 15,16 combined with an accecptance-rejection method 14 we implement and demonstrate the WOH algorithm for the induced charge density distribution on parallel infinite conducting plates when a unit charge is between the plates and apply it to the mutual capacitance of two circular parallel plates. In both simulations, WOH algorithm shows much better performance than the WOS algorithm.
Diffusion Monte Carlo algorithms have been mainly applied for extracting mutual capacitances for a system of conductors, 12,13,31,32 . In semiconductor industry, a powerful commercial 3D CAD tool, QuickCap TM has been used 33 . It is emphasized that WOH algorithm can be used for the cases of diffusion near to a finite plane boundary or the one surrounded by two-plane between boundaries.

Data availability
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.  Table 2. Runtime (in seconds) comparison of WOH algorithm with WOS one; for the mutual capacitance of the two parallel circular plates we use 100 independent simulations of 10 10 MC steps. R/D is the ratio of the radius of the circular plate and the separation between the plates.